On the danger of detecting network states in white noise

The general idea of nonstationarity of brain activity or dependence of the dynamics on some, potentially unobserved, temporally changing or fluctuating parameter, has been familiar in the neuroscience community in contexts such as sleep dynamics or epileptology for a long time. However, recently it has been attracting increasing attention in the context of functional brain network analysis. This seems as a natural development of the field—once that functional connectivity as computed under the simplifying stationarity assumption has been well established, it is only logical to try to detect changes in brain functional connectivity over time. In general, detecting such nonstationarities in a reliable fashion is a methodologically challenging task, as changes in estimates of functional connectivity over time may be also due to random fluctuations, rather than genuine changes of the process. There is a wide array of approaches to studying such nonstationarities documented in literature (Hutchison et al., 2013), and an important but often neglected general methodological step is assessing the results against an appropriate null model corresponding to stationary process. In the following, we give an illustrative example of how a typical nonstationarity analysis can generate spurious signs of nonstationary dynamics even when applied to stationary process. To show that this is not a purely theoretical issue, we closely follow the analysis procedure used in a recently published study by Betzel et al. (2012). We note that this particular paper have caught our attention by coincidence, while we believe the issue is pertinent to a substantial fraction of the literature. In their paper, (Betzel et al., 2012) deal with characterizing the dynamics of brain activity measured by EEG. In particular, Betzel et al. report the detection of rapid transitions between intermittently stable states, explicitly saying that “As predicted, fast (∼100 ms) dynamics of wholebrain synchronization were observed during resting-state EEG,” documenting the typical fast (∼100ms) time scale of these states in Figure 6B of their paper (see also Figures 4, 5). Their argument is based on the following data-processing scheme: First, for each time point of filtered EEG data, a functional connectivity matrix is computed using pairwise synchronization likelihood values and the time points are clustered based on similarity of the corresponding functional connectivity matrices. Next, contiguous stretches of time points that are members of the same cluster are interpreted as corresponding to a duration of an atomic brain state. Finally, the brain-state-representing functional connectivity matrices are pooled across subjects and clustered based on their similarity to define higher-order states. Notably, the procedure applied by Betzel et al. is principally data-driven, rather than relying on some model testing or assumptions, and it includes band-pass filtering and sliding-window-like analysis. We therefore conjectured that the temporal structure of the observed functional connectivity dynamics might have been crucially affected by the procedure itself (as the authors tentatively admitted in their discussion, albeit unfortunately have not tested the results against stationary model data). To explore the viability of this alternative explanation, we applied a processing pipeline built according to the description given in the original manuscript to model data, consisting of 100 samples (each of length T = 2500 time points, representing mock 5 s epoch of EEG data) of a multivariate (N = 20) white noise process. The applied processing steps included application of frequency filtering (using elliptic filters corresponding to the four specified frequency bands; we applied zero-phase digital filtering by processing the input data in both the forward and reverse directions) and subsequent computation of the synchronization likelihood (Stam and van Dijk, 2002). The parameters of the synchronization likelihood l, m, w1, w2, nrec were set for each frequency band as in Betzel et al. (2012). The resulting functional connectivity matrices were clustered using the standard k-means clustering method (Lloyd, 1982). In Figure 1 you see that the typical duration of detected states closely corresponds to the distributions observed in the original paper (compare with Figures 6B, 4A,B in Betzel et al., 2012). In particular, the typical timescale is in the order of tens to hundreds of ms. Also, this time scale depends on the selected filtering in the same way as in the original work, with the time scales of the beta and theta bands